from cmath import cos
from lib2to3.pgen2.pgen import DFAState
import numpy as np
import pandas as pd
import os
import sys
from matplotlib import pyplot as plt


def GUM_RayD(data_L2,
             data_L0,
             param='D',
             z0=30,
             dens0=0.01841,
             u_dens0=0,
             tau=1e-9,
             u_tau=0,
             alt_angle=0,
             location_altitude=0,
             height_resolution=1,
             time_resolution=3600,
             pulse_repeat=30,
             k=1):
    # data_L0：L0级数据，第1列为高度，第2列为光子数
    # data_L2：L2级数据，第1列为高度，第2列为密度
    # tau：探测器死区时间
    # u_tau：探测器死区时间波动度
    # pulse_repeat：激光器重频
    # time_resolution：L2级积分时间
    # location_altitude：系统所在高度
    # z0：参考密度所在高度
    # dens0：参考密度
    # uncertain_densityRef：参考密度波动度
    # tau = 1e-9
    # u_tau = 1e-9 * 0.01
    # time_resolution = 1200
    # pulse_repeat = 30
    # z0 = 30  #参考密度所在高度
    # dens0 = 3.64686762151e+017  #参考密度
    # u_dens0 = 0  #参考密度波动度
    # location_altitude = 1.14
    z = data_L2[:, 0]
    dens = data_L2[:, 1]
    #海拔高度转化和距离门内光子数累加
    data_L0_s = get_bins_data(data_L0, z, location_altitude, alt_angle, height_resolution)
    #背景噪声获取
    B, B_mean_num = get_backgroud_noise(data_L0_s)
    B = B[0]
    #探测范围截取
    Nr1 = get_target_data(data_L0_s, z)[:, 0]  #目标高度光子数
    Nr2 = get_target_data(data_L0_s, z0)[0]  #参考高度处光子数
    # z = z[:np.size(Nr1)]
    # dens = dens[:np.size(Nr1)]
    #光子数校正过程中的灵敏系数
    Nc1 = Nr1
    Nb1 = Nc1 - B
    allSize = time_resolution * pulse_repeat  #L2数据累加的脉冲数
    c_Nb1_tau1 = Nc1**2 / allSize
    c_Nb1_Nr1 = Nc1**2 / Nr1**2
    c_Nb1_B = 1
    Nc2 = Nr2
    Nb2 = Nc2 - B
    c_Nb2_tau2 = Nc2**2 / allSize
    c_Nb2_Nr2 = Nc2**2 / Nr2**2
    c_Nb2_B = 1
    #反演过程中灵敏系数
    c_dens_Nb1 = dens / Nb1
    c_dens_Nb2 = dens / Nb2
    c_dens_dens0 = dens / dens0
    #光子噪声不确定度分量
    u_Nr1 = np.sqrt(Nr1)
    u_Nr2 = np.sqrt(Nr2)
    c_Nr1 = c_dens_Nb1 * c_Nb1_Nr1
    c_Nr2 = c_dens_Nb2 * c_Nb2_Nr2
    u1 = np.abs(u_Nr1 * c_Nr1)
    u2 = np.abs(u_Nr2 * c_Nr2)
    #背景噪声不确定度分量
    u_B1 = np.sqrt(B / B_mean_num)
    c_B1 = c_dens_Nb1 * c_Nb1_B
    c_B2 = c_dens_Nb2 * c_Nb2_B
    c_B = c_B1 + c_B2
    u3 = np.abs(u_B1 * c_B)
    #死区时间不确定度分量
    c_tau1 = c_dens_Nb1 * c_Nb1_tau1
    c_tau2 = c_dens_Nb2 * c_Nb2_tau2
    u4 = np.abs(u_tau * c_tau1)
    u5 = np.abs(u_tau * c_tau2)
    #参考密度不确定度分量
    u6 = np.abs(u_dens0 * c_dens_dens0)
    #合成标准不确定度
    u = np.sqrt(u1**2 + u2**2 + u3**2 + u4**2 + u5**2 + u6**2)
    #扩展不确定度
    U = k * u
    #不确定度相对百分比
    U = U / dens * 100
    result = np.hstack((z.reshape(-1, 1), dens.reshape(-1, 1), U.reshape(-1, 1)))
    if param == 'D':
        unit = 'kg/m^3'
    else:
        unit = '1/cm^3'
    header = '高度(km) 密度({}) U(%)_MCM标准不确定度'.format(
        unit)
    result_df = pd.DataFrame(result, columns=header.split(' '))
    # result_df.drop('U(%)_GUM扩展不确定度', axis=1, inplace=True)
    return result_df


def GUM_RayT(data_L2,
             data_L0,
            z0=30,
             dens0=0.01841,
             u_dens0=0,
             z1=80,
             T1=209.6083,
             u_T1=0,
             tau=1e-9,
             u_tau=0,
             alt_angle=0,
             location_altitude=0,
             height_resolution=1,
             time_resolution=3600,
             pulse_repeat=30,
             k=1):
    # data_L0：L0级数据，第1列为高度，第2列为RayVHS
    # data_L2：L0级数据，第1列为高度，第2列为RayVHS
    # tau：探测器死区时间
    # u_tau：探测器死区时间波动度
    # pulse_repeat：激光器重频
    # time_resolution：L2级积分时间
    # location_altitude：系统所在高度
    # z0：参考密度所在高度
    # dens0：参考密度
    # uncertain_densityRef：参考密度波动度
    # tau = 1e-9
    # u_tau = 1e-9 * 0.01
    # time_resolution = 1200
    # pulse_repeat = 30
    # z0 = 30  #参考密度所在高度
    # dens0 = 3.64686762151e+017  #参考密度
    # u_dens0 = 0  #参考密度波动度
    # z1 = 80  #密度积分起始高度
    # T1 = 209.6083 #密度积分起始温度
    # u_T1 = 0  #参考温度波动度
    z = data_L2[:, 0]
    T = data_L2[:, 1]
    # todo
    # mask = z <= z1
    # z = z[mask]
    # T = T[mask]
    #海拔高度转化和距离门内光子数累加
    data_L0_s = get_bins_data(data_L0, z, location_altitude, alt_angle, height_resolution)
    #背景噪声获取
    B, B_mean_num = get_backgroud_noise(data_L0_s)
    B = B[0]
    #探测范围截取
    Nr = get_target_data(data_L0_s, z)[:, 0]  #目标高度光子数
    Nr0 = get_target_data(data_L0_s, z0)[0]  #参考高度处光子数
    Nr1 = get_target_data(data_L0_s, z1)[0]  #参考高度处光子数
    # np.savetxt('L0.txt', data_L0_s)
    #光子数校正过程中的灵敏系数
    Nc = Nr
    Nb = Nc - B
    allSize = time_resolution * pulse_repeat  #L2数据累加的脉冲数
    c_Nb_tau1 = Nc**2 / allSize
    c_Nb_Nr = Nc**2 / Nr**2
    c_Nb_B = 1
    Nc0 = Nr0
    Nb0 = Nc0 - B
    c_Nb0_tau2 = Nc0**2 / allSize
    c_Nb0_Nr0 = Nc0**2 / Nr0**2
    c_Nb0_B = 1
    Nc1 = Nr1
    Nb1 = Nc1 - B
    c_Nb1_tau2 = Nc1**2 / allSize
    c_Nb1_Nr1 = Nc1**2 / Nr1**2
    c_Nb1_B = 1
    #反演过程中灵敏系数
    c_T_T1 = Nb1 / Nb * z1**2 / z**2
    c_T_Nb = -T / Nb
    c_T_Nb0 = T / Nb0
    #光子噪声不确定度分量
    u_Nr = np.sqrt(Nr)
    u_Nr0 = np.sqrt(Nr0)
    c_Nr = c_T_Nb * c_Nb_Nr
    c_Nr0 = c_T_Nb0 * c_Nb0_Nr0
    u1 = np.abs(u_Nr * c_Nr)
    u2 = np.abs(u_Nr0 * c_Nr0)
    #背景噪声不确定度分量
    u_B = np.sqrt(B / B_mean_num)
    c_B1 = c_T_Nb * c_Nb_B
    c_B2 = c_T_Nb0 * c_Nb0_B
    c_B = c_B1 + c_B2
    u3 = np.abs(u_B * c_B)
    #死区时间不确定度分量
    c_tau1 = c_T_Nb * c_Nb_tau1
    c_tau2 = c_T_Nb0 * c_Nb0_tau2
    u4 = np.abs(u_tau * c_tau1)
    u5 = np.abs(u_tau * c_tau2)
    #参考温度不确定度分量
    u6 = np.abs(u_T1 * c_T_T1)
    #合成标准不确定度
    u = np.sqrt(u1**2 + u2**2 + u3**2 + u4**2 + u5**2 + u6**2)
    #扩展不确定度
    U = k * u
    #不确定度相对百分比
    # U, u, u1, u2, u3, u4, u5, u6, u7, u8 = (U, u, u1, u2, u3, u4, u5, u6, u7, u8) / dens * 100
    result = np.hstack((z.reshape(-1, 1), T.reshape(-1, 1), U.reshape(-1, 1), u.reshape(-1, 1), u1.reshape(-1, 1),
                        u2.reshape(-1, 1), u3.reshape(-1, 1), u4.reshape(-1, 1), u5.reshape(-1, 1), u6.reshape(-1, 1)))
    header = '高度(km) 温度(K) U(K)_GUM扩展不确定度 u(K)_标准不确定度 u1(K)_探测点光子噪声 u2(K)_参考点光子噪声 u3(K)_背景噪声 u4(K)_探测点死区时间 u5(K)_参考点死区时间 u6(K)_参考温度'
    result_df = pd.DataFrame(result, columns=header.split(' '))
    # result_df.drop('U(K)_GUM扩展不确定度', axis=1, inplace=True)
    return result_df


def GUM_RayW(data_L2,
             data_L0,
             mode=0,
             c_V_RV=707,
             wavelen=532,
             u_wavelen=0.01,
             tau=1e-9,
             u_tau=0,
             alt_angle=0,
             location_altitude=0,
             height_resolution=1,
             time_resolution=3600,
             pulse_repeat=30,
             k=1):
    # c_V_RV：测风灵敏度
    # mode: 0单边缘，1双边缘

    # data_L0：L0级数据，第1列为高度，第2列为光子数
    # data_L2：L2级数据，第1列为高度，第2列为密度
    # tau：探测器死区时间
    # u_tau：探测器死区时间波动度
    # pulse_repeat：激光器重频
    # time_resolution：L2级积分时间
    # location_altitude：系统所在高度
    # z0：参考密度所在高度
    # dens0：参考密度
    # uncertain_densityRef：参考密度波动度
    # tau = 1e-9
    # u_tau = 1e-9 * 0.01
    # time_resolution = 1200
    # pulse_repeat = 30
    # z0 = 30  #参考密度所在高度
    # dens0 = 3.64686762151e+017  #参考密度
    # u_dens0 = 0  #参考密度波动度
    # location_altitude = 1.14
    z = data_L2[:, 0]
    V = data_L2[:, 1]
    #海拔高度转化和距离门内光子数累加
    data_L0_s = get_bins_data(data_L0, z, location_altitude, alt_angle, height_resolution)
    #背景噪声获取
    B, B_mean_num = get_backgroud_noise(data_L0_s)
    B = B[0]
    #探测范围截取
    Nr = get_target_data(data_L0_s, z)  #目标高度光子数
    Nr1 = Nr[:, 0]
    Nr2 = Nr[:, 1]
    #光子数校正过程中的灵敏系数
    Nc1 = Nr1
    Nb1 = Nc1 - B
    allSize = time_resolution * pulse_repeat  #L2数据累加的脉冲数
    c_Nb1_tau1 = Nc1**2 / allSize
    c_Nb1_Nr1 = Nc1**2 / Nr1**2
    c_Nb1_B = 1
    Nc2 = Nr2
    Nb2 = Nc2 - B
    c_Nb2_tau2 = Nc2**2 / allSize
    c_Nb2_Nr2 = Nc2**2 / Nr2**2
    c_Nb2_B = 1
    #反演过程中灵敏系数
    if mode == 0:
        # 针对RV = (N1-N2)/(N1+N2)的情況
        # if c_V_RV == None:
        #     RV = (Nb1 - Nb2) / (Nb1 + Nb2)
        #     c_V_RV = np.polyfit(RV, V, 1)[0]
        c_RV_Nb1 = 2 * Nb2 / (Nb1 + Nb2)**2
        c_RV_Nb2 = 2 * Nb1 / (Nb1 + Nb2)**2
    else:
        # 针对RV = N1/N2的情況
        # if c_V_RV == None:
        #     RV = Nb1 / Nb2
        #     c_V_RV = np.polyfit(RV, V, 1)[0]
        c_RV_Nb1 = 1 / Nb2
        c_RV_Nb2 = Nb1 / Nb2**2
    # plt.figure()
    # plt.scatter(RV, V)
    # plt.plot(RV, RV * c_V_RV + np.polyfit(RV, V, 1)[1])
    # plt.savefig('ResultPic/风速偏导拟合.png')
    # plt.close()
    c_V_Nb1 = c_V_RV * c_RV_Nb1
    c_V_Nb2 = c_V_RV * c_RV_Nb2
    # c_V_V0 = V / V0
    #光子噪声不确定度分量
    u_Nr1 = np.sqrt(Nr1)
    u_Nr2 = np.sqrt(Nr2)
    c_Nr1 = c_V_Nb1 * c_Nb1_Nr1
    c_Nr2 = c_V_Nb2 * c_Nb2_Nr2
    u1 = np.abs(u_Nr1 * c_Nr1)
    u2 = np.abs(u_Nr2 * c_Nr2)
    #背景噪声不确定度分量
    u_B1 = np.sqrt(B / B_mean_num)
    c_B1 = c_V_Nb1 * c_Nb1_B
    c_B2 = c_V_Nb2 * c_Nb2_B
    c_B = c_B1 + c_B2
    u3 = np.abs(u_B1 * c_B)
    #死区时间不确定度分量
    c_tau1 = c_V_Nb1 * c_Nb1_tau1
    c_tau2 = c_V_Nb2 * c_Nb2_tau2
    u4 = np.abs(u_tau * c_tau1)
    u5 = np.abs(u_tau * c_tau2)
    #参考波长不确定度分量
    c_wavelen = 3e8 / (wavelen * 1e-9) / 2
    u6 = np.abs(u_wavelen * c_wavelen) * np.ones(np.shape(V))
    #合成标准不确定度
    u = np.sqrt(u1**2 + u2**2 + u3**2 + u4**2 + u5**2 + u6**2)
    #扩展不确定度
    U = k * u
    result = np.hstack((z.reshape(-1, 1), V.reshape(-1, 1), U.reshape(-1, 1), u.reshape(-1, 1), u1.reshape(-1, 1),
                        u2.reshape(-1, 1), u3.reshape(-1, 1), u4.reshape(-1, 1), u5.reshape(-1, 1), u6.reshape(-1, 1)))
    header = '高度(km) 风速(m/s) U(m/s)_GUM扩展不确定度 u(m/s)_标准不确定度 u1(m/s)_探测通道光子噪声 u2(m/s)_参考通道光子噪声 u3(m/s)_背景噪声 u4(m/s)_探测通道死区时间 u5(m/s)_参考通道死区时间 u6(m/s)_参考波长'
    result_df = pd.DataFrame(result, columns=header.split(' '))
    # result_df.drop('U(%)_GUM扩展不确定度', axis=1, inplace=True)
    return result_df


def GUM_ResT(data_L2,
             data_L0,
             mode=0,
             c_T_RT=0.00349,
             c_T_RV=0.00349,
             tau=1e-9,
             u_tau=0,
             alt_angle=0,
             location_altitude=0,
             height_resolution=1,
             time_resolution=3600,
             pulse_repeat=30,
             k=1):
    # c_T_RT: 测温灵敏度，dT/dR_T
    # c_V_RV：测风灵敏度，dV/dR_V
    # data_L0：L0级数据，第1列为高度，第2列为RayVHS
    # data_L2：L2级数据，第1列为高度，第2列为密度TempRay
    # tau：探测器死区时间
    # u_tau：探测器死区时间波动度
    # pulse_repeat：激光器重频
    # time_resolution：L2级积分时间
    # location_altitude：系统所在高度
    # z0：参考密度所在高度
    # dens0：参考密度
    # uncertain_densityRef：参考密度波动度
    # tau = 1e-9
    # u_tau = 1e-9 * 0.01
    # time_resolution = 1200
    # pulse_repeat = 30
    # z0 = 30  #参考密度所在高度
    # dens0 = 3.64686762151e+017  #参考密度
    # u_dens0 = 0  #参考密度波动度
    # location_altitude = 1.14
    z = data_L2[:, 0]
    T = data_L2[:, 1]
    #海拔高度转化和距离门内光子数累加
    data_L0_s = get_bins_data(data_L0, z, location_altitude, alt_angle, height_resolution)
    #背景噪声获取
    B, B_mean_num = get_backgroud_noise(data_L0_s)
    B1 = B[0]
    B2 = B[1]
    B3 = B[2]
    #探测范围截取
    Nr = get_target_data(data_L0_s, z)  #目标高度光子数
    Nr1 = Nr[:, 0]
    Nr2 = Nr[:, 1]
    Nr3 = Nr[:, 2]
    #光子数校正过程中的灵敏系数
    Nc1 = Nr1
    Nb1 = Nc1 - B1
    allSize = time_resolution * pulse_repeat  #L2数据累加的脉冲数
    c_Nb1_tau1 = Nc1**2 / allSize
    c_Nb1_Nr1 = Nc1**2 / Nr1**2
    c_Nb1_B = 1
    Nc2 = Nr2
    Nb2 = Nc2 - B2
    c_Nb2_tau2 = Nc2**2 / allSize
    c_Nb2_Nr2 = Nc2**2 / Nr2**2
    c_Nb2_B = 1
    Nc3 = Nr3
    Nb3 = Nc3 - B3
    c_Nb3_tau3 = Nc3**3 / allSize
    c_Nb3_Nr3 = Nc3**3 / Nr3**3
    c_Nb3_B = 1
    #反演过程中灵敏系数
    if mode == 0:
        # 针对RT=(Nb1+Nb2)/Nb3/2, RV=(Nb1-Nb2)/Nb3
        c_RT_Nb1 = 1 / Nb3 / 2
        c_RV_Nb1 = 1 / Nb3
        c_RT_Nb2 = 1 / Nb3 / 2
        c_RV_Nb2 = -1 / Nb3
        c_RT_Nb3 = -(Nb1 + Nb2) / Nb3**2 / 2
        c_RV_Nb3 = -(Nb1 - Nb2) / Nb3**2
        # if c_T_RT == None:
        #     RT = (Nb1 + Nb2) / Nb3 / 2
        #     c_T_RT = np.polyfit(RT, T, 1)[0]
        # if c_T_RV == None:
        #     RV = (Nb1 - Nb2) / Nb3
        #     c_T_RV = np.polyfit(RV, T, 1)[0]
    c_T_Nb1 = c_T_RT * c_RT_Nb1 + c_T_RV * c_RV_Nb1
    c_T_Nb2 = c_T_RT * c_RT_Nb2 + c_T_RV * c_RV_Nb2
    c_T_Nb3 = c_T_RT * c_RT_Nb3 + c_T_RV * c_RV_Nb3
    #光子噪声不确定度分量
    u_Nr1 = np.sqrt(Nr1)
    u_Nr2 = np.sqrt(Nr2)
    u_Nr3 = np.sqrt(Nr3)
    c_Nr1 = c_T_Nb1 * c_Nb1_Nr1
    c_Nr2 = c_T_Nb2 * c_Nb2_Nr2
    c_Nr3 = c_T_Nb3 * c_Nb3_Nr3
    u1 = np.abs(u_Nr1 * c_Nr1)
    u2 = np.abs(u_Nr2 * c_Nr2)
    u3 = np.abs(u_Nr3 * c_Nr3)
    #背景噪声不确定度分量
    u_B1 = np.sqrt(B1 / B_mean_num)
    u_B2 = np.sqrt(B2 / B_mean_num)
    u_B3 = np.sqrt(B3 / B_mean_num)
    c_B1 = c_T_Nb1 * c_Nb1_B
    c_B2 = c_T_Nb2 * c_Nb2_B
    c_B3 = c_T_Nb3 * c_Nb3_B
    u4 = np.abs(u_B1 * c_B1)
    u5 = np.abs(u_B2 * c_B2)
    u6 = np.abs(u_B3 * c_B3)
    #死区时间不确定度分量
    c_tau1 = c_T_Nb1 * c_Nb1_tau1
    c_tau2 = c_T_Nb2 * c_Nb2_tau2
    c_tau3 = c_T_Nb2 * c_Nb3_tau3
    u7 = np.abs(u_tau * c_tau1)
    u8 = np.abs(u_tau * c_tau2)
    u9 = np.abs(u_tau * c_tau3)
    #合成标准不确定度
    u = np.sqrt(u1**2 + u2**2 + u3**2 + u4**2 + u5**2 + u6**2 + u7**2 + u8**2 + u9**2)
    #扩展不确定度
    U = k * u
    # #不确定度相对百分比
    # U, u, u1, u2, u3, u4, u5, u6 = (U, u, u1, u2, u3, u4, u5, u6) / dens * 100
    result = np.vstack((z, T, U, u, u1, u2, u3, u4, u5, u6, u7, u8, u9)).T
    header = '高度(km) 温度(K) U(K)_GUM扩展不确定度 u(K)_标准不确定度 u1(K)_峰值频率通道光子噪声 u2(K)_右翼频率通道光子噪声 u3(K)_左翼频率通道光子噪声 u4(K)_峰值频率通道背景噪声 u5(K)_右翼频率通道背景噪声 u6(K)_左翼频率通道背景噪声 u7(K)_峰值频率通道死区时间 u8(K)_右翼频率通道死区时间 u9(K)_左翼频率通道死区时间'
    result_df = pd.DataFrame(result, columns=header.split(' '))
    return result_df


def GUM_ResW(data_L2,
             data_L0,
             mode=0,
             c_V_RT=0.00349,
             c_V_RV=0.00349,
             tau=1e-9,
             u_tau=0,
             alt_angle=0,
             location_altitude=0,
             height_resolution=1,
             time_resolution=3600,
             pulse_repeat=30,
             k=1):
    pass


def GUM_RamanT(data_L2,
               data_L0,
               poly_coef,
               mode=0,
               tau=1e-9,
               u_tau=0,
               alt_angle=0,
               location_altitude=0,
               height_resolution=1,
               time_resolution=3600,
               pulse_repeat=30,
               k=1):
    #输入量
    poly_coef = [0.003945, -0.00171, -0.00101]
    a = poly_coef[0]
    b = poly_coef[1]
    c = poly_coef[2]
    z = data_L2[:, 0]
    T = data_L2[:, 1]
    #海拔高度转化和距离门内光子数累加
    data_L0_s = np.copy(data_L0)
    data_L0_s[:, 0] = data_L0_s[:, 0] * np.cos(alt_angle / 180 * np.pi) + location_altitude
    dz = np.abs(np.diff(z).mean())
    data_L0_s[:, 0] = np.round((data_L0_s[:, 0] - z[0]) / dz) * dz + z[0]
    data_L0_s = pd.DataFrame(data_L0_s).groupby(0, as_index=False).sum().values
    #背景噪声获取
    B1 = data_L0_s[-1, 1]
    B2 = data_L0_s[-1, 2]
    #探测范围截取
    ind = np.nanargmin(np.abs(z - z0))
    Nr2 = data_L0_s[ind, 1]  #参考高度处光子数
    mask = (z[0] <= data_L0_s[:, 0]) & (data_L0_s[:, 0] <= z[-1])
    data_L0_s = data_L0_s[mask, :]
    Nr1 = data_L0_s[:, 1]
    #光子数校正过程中的灵敏系数
    Nc1 = Nr1
    Nb1 = Nc1 - B1
    allSize = time_resolution * pulse_repeat  #L2数据累加的脉冲数
    c_Nb1_tau1 = Nc1**2 / allSize
    c_Nb1_Nr1 = Nc1**2 / Nr1**2
    c_Nb1_B1 = 1
    Nc2 = Nr2
    Nb2 = Nc2 - B2
    c_Nb2_tau2 = Nc2**2 / allSize
    c_Nb2_Nr2 = Nc2**2 / Nr2**2
    c_Nb2_B2 = 1
    #反演过程中灵敏系数
    Q = Nb1 / Nb2
    temp = 1 / (c * np.log(Q)**2 + b * np.log(Q) + a)
    c_T_Nb1 = 1 / Nc2 * temp**2 * (2 * c * np.log(Q) + b) / Q
    c_T_Nb2 = Nc1 / Nc2**2 * temp**2 * (2 * c * np.log(Q) + b) / Q
    #光子噪声不确定度分量
    u_Nr1 = np.sqrt(Nr1)
    u_Nr2 = np.sqrt(Nr2)
    c_Nr1 = c_T_Nb1 * c_Nb1_Nr1
    c_Nr2 = c_T_Nb2 * c_Nb2_Nr2
    u1 = np.abs(u_Nr1 * c_Nr1)
    u2 = np.abs(u_Nr2 * c_Nr2)
    #背景噪声不确定度分量
    u_B1 = np.sqrt(B1)
    u_B2 = np.sqrt(B2)
    c_B1 = c_T_Nb1 * c_Nb1_B1
    c_B2 = c_T_Nb2 * c_Nb2_B2
    u3 = np.abs(u_B1 * c_B1)
    u4 = np.abs(u_B2 * c_B2)
    #死区时间不确定度分量
    c_tau1 = c_T_Nb1 * c_Nb1_tau1
    c_tau2 = c_T_Nb2 * c_Nb2_tau2
    u5 = np.abs(u_tau * c_tau1)
    u6 = np.abs(u_tau * c_tau2)
    #合成标准不确定度
    u = np.sqrt(u1**2 + u2**2 + u3**2 + u4**2 + u5**2 + u6**2)


def get_bins_data(data_L0, z, location_altitude=0, alt_angle=0, height_resolution=1):
    z = np.round(z, 3)
    data_L0_s = np.copy(data_L0)
    data_L0_s[:, 0] = data_L0_s[:, 0] * np.sin(alt_angle / 180 * np.pi) + location_altitude
    dz = np.abs(np.diff(z).mean())
    data_L0_s[:, 0] = np.round((data_L0_s[:, 0] - z[0]) / dz) * dz + z[0]
    data_L0_s[:, 0] = np.round(data_L0_s[:, 0], 3)
    data_L0_s = pd.DataFrame(data_L0_s).groupby(0, as_index=False).sum().values
    if height_resolution == None:
        data_L0_s[:, 1:] = data_L0_s[:, 1:] * height_resolution / dz
    return data_L0_s


def get_backgroud_noise(data_L0_s):
    mask = (140 <= data_L0_s[:, 0]) & (data_L0_s[:, 0] <= 160)
    B_mean_num = np.size(mask)
    B = np.nanmean(data_L0_s[mask, 1:], axis=0)  # 背景噪声取法1,以140km到160km均值为背景噪声
    # B = data_L0_s[-1, 1]  # 背景噪声取法2,以最后一个距离门为背景噪声
    return B, B_mean_num


def get_target_data(data_L0_s, z):
    z = np.round(z, 3)
    if isinstance(z, (list, tuple, np.ndarray)):
        mask = (z[0] <= data_L0_s[:, 0]) & (data_L0_s[:, 0] <= z[-1])
        data_L0_s = data_L0_s[mask, :]
        Nr = data_L0_s[:, 1:]
    else:
        ind = np.nanargmin(np.abs(data_L0_s[:, 0] - z))
        Nr = data_L0_s[ind, 1:]  #参考高度z处光子数
    return Nr
